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ABSTRACT 


Recent studies have revealed a strong relation between sample-averaged black-hole (BH) accretion rate 
(BHAR) and star formation rate (SFR) among bulge-dominated galaxies, i.e., "lockstep" BH-bulge growth, 
in the distant universe. This relation might be closely related to the BH-bulge mass correlation observed in the 
local universe. To understand further BH-bulge coevolution, we present ALMA CO(2-1) or CO(3-2) observa- 
tions of 7 star-forming bulge-dominated galaxies at z = 0.5-2.5. Using the ALMA data, we detect significant 
(> 3c) CO emission from 4 objects. For our sample of 7 galaxies, we measure (or constrain with upper limits) 
their CO line fluxes and estimate molecular gas masses (Mgas). We also estimate their stellar masses (Mstar) 
and SFRs by modelling their spectral energy distributions (SEDs). Using these physical properties, we derive 
the gas-depletion timescales (rae = Mgas/SFR) and compare them with the bulge/BH growth timescales 
(Terow = Mstar /SFR ~ Mgu/BHAR). Our sample generally has Tdep Shorter than Tgrow by a median factor 
of = 4, indicating that the cold gas will be depleted before significant bulge/BH growth takes place. This re- 
sult suggests that the BH-bulge lockstep growth is mainly responsible for maintaining their mass relation, not 
creating it. We note that our sample is small and limited to z « 2.5; JWST and ALMA will be able to probe to 


higher redshifts in the near future. 
1. INTRODUCTION 


From observations of the nearby universe, massive galax- 
ies always host supermassive black holes (BHs) in their cen- 
tral regions, and the BH mass (Mgg) is tightly correlated 
with the stellar mass (Mstar) of the host-galaxy bulge (e.g., 
Kormendy & Ho 2013; Saglia et al. 2016). This tight BH- 
bulge mass correlation indicates that a physical connection 
between the BHs and bulges exists in their cosmic evolution 
history. This connection is often called “BH-bulge coevolu- 
tion". 

The specific coevolution mechanisms remain largely un- 
known. Some early studies proposed minor-merger events 
as the driver of the tight BH-bulge mass relation in the lo- 
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cal universe (e.g., Peng 2007; Jahnke & Macció 2011). The 
idea is based on the statistical central-limit theorem: if low- 
mass galaxies are scattered around the Mpuy-Mstar relation, 
the final system will be close to the mass relation after many 
episodes of minor mergers. This scheme requires that low- 
mass systems are centered around the Mpy-Msgtar relation. 
However, more recent observations indicate that low-mass 
galaxies have Mgg systematically below the relation (e.g., 
Aird et al. 2018; Yang et al. 2018) and that some low-mass 
galaxies may not even host BHs (e.g., Greene et al. 2020). 
Therefore, minor mergers are unlikely to be the origin of the 
BH-bulge mass relation. Another idea about the origin of 
the Mpy-Mgtar relation is related to AGN feedback. A BH 
grows to a critical mass (related to Mpuige) and launches a 
powerful wind that removes cold gas and/or prevents gas re- 
plenishment (e.g., King & Pounds 2015). The BH and stel- 
lar growth are thereby halted due to the lack of fuel. This 


2 YANG ET AL. 


scenario requires strong AGN negative feedback on star for- 
mation. However, this assumption still lacks strong observa- 
tional support (e.g., Harrison 2017; Shangguan et al. 2020). 
From an observational point of view, the Mpy-Mstar re- 
lation might result from some connections between BH ac- 
cretion and host galaxies. Earlier observations on this topic 
mostly focused on the connections of BH growth versus 
host-galaxy Mstar and star formation rate (SFR) at differ- 
ent redshifts (see $3.3 of Brandt & Yang 2021 for a review). 
Some significant relations have been found. For example, 
BH growth strongly depends upon Mstar at a given redshift 
(e.g., Georgakakis et al. 2017; Yang et al. 2017, 2018; Aird 
et al. 2018). However, it is not obvious how these relations 
are related to the BH-bulge correlation in the local universe, 
mainly because these relations are for global galaxy (rather 
than bulge) properties. It is technically challenging to inves- 
tigate bulges, since galaxy angular sizes are generally small 
(S arcsec scales) in the distant universe (z = 1). Also, to 
avoid strong effects from the morphological K correction" 
in the rest-frame UV, the sampled rest-frame wavelengths are 
required to be longer than ~ 4000A (e.g., Papovich et al. 
2005; Conselice 2014). This means that the imaging should 
use IR bands (= 1.2 um) to measure the morphologies of 
sources at z ~ 2, the cosmic epoch when active galactic nu- 
cleus (AGN) and star formation (SF) activities peak. The 
Cosmic Assembly Near-infrared Deep Extragalactic Legacy 
Survey (CANDELS; Grogin et al. 2011; Koekemoer et al. 
2011) provides high-angular-resolution (~ 0.2") and deep 
H-band imaging over a ~ 900 arcmin? area, providing an 
excellent chance to study AGN-bulge connections at z < 3. 
Using the CANDELS-based morphological measurements 
(Huertas-Company et al. 2015) as well as other multiwave- 
length data, Yang et al. (2019) selected a sample of bulge- 
dominated galaxies and investigated their sample-averaged 
BH accretion rate (BHAR) versus SFR.! They found a sig- 
nificant linear correlation between BHAR and SFR for the 
bulge-dominated galaxies, and this correlation does not ex- 
ist among their comparison sample that has disks and/or ir- 
regularities (see, e.g., Kocevski et al. 2017; Ni et al. 2019, 
2021 for related results). The best-fit BHAR/SFR ratio in 
Yang et al. (2019) is ~ 1/300, similar to the BH/bulge mass 
ratio observed in the local universe (e.g., Kormendy & Ho 
2013). This “lockstep” style growth between BHs and bulges 
can be useful for predicting BH accretion from bulge SF in- 
formation. Based on a sample of bulge-dominated galaxies 
with well-measured star formation histories (SFHs; Estrada- 


This sample-averaged BHAR is designed to overcome AGN short-term 
(€ 107 years) variability effects and approximate long-term BH growth 
rate (e.g., Hickox et al. 2014). Therefore, individual bulge-dominated sys- 
tems should follow the same BHAR-SFR relation over cosmic evolution 
timescales (> 108 years). 


Carpenter et al. 2020), Yang et al. (2021) predicted BH ac- 
cretion densities at high redshifts, which can be tested with 
JWST and future IR missions. 

The BHAR-bulge SFR relation could be closely related to 
the BH-bulge mass relation in the local universe, given their 
similarities as above. However, the detailed picture of BH- 
bulge coevolution is still not clear, and there are two possible 
scenarios following the formation of the bulge: 


1. The system initially retains a large amount of cold 
gas, which is available to fuel significant bulge growth. 
The BH then would continue to grow following the 
BHAR-SFR relation. The galaxy thereby moves to 
Mpu/Mstar ~ BHAR/SFR z 1/300, regardless 
of the previous ratio of Mpy/Mgtar in the pre-bulge 
phase. In this case, the BHAR-SFR relation plays a 
role of "creation" for the BH-bulge mass correlation. 


2. The system initially has a limited or insignificant 
amount of cold gas remaining. Star formation in the 
bulge will soon shut down soon due to the lack of fuel, 
assuming no further gas replenishment. The galaxy 
thereby remains in a low-specific SFR state and has 
a low-BHAR (according to the BHAR-SFR relation), 
and it passively evolves. The Mpy/Msgtar “freezes” 
at the value in the pre-bulge phase. In this case, the 
BHAR-SFR relation plays a role of “maintenance” for 
the BH-bulge mass correlation. 


Fig. | illustrates these two scenarios above schematically. 
The determining factor between the two scenarios is the 
amount of cold gas available after bulge formation. If gas 
is abundant, then the more realistic scenario is the first one. 
The system has significant BH/bulge mass growth, reducing 
the scatter around the local Mpy-Mstar relation (Fig. 1 left). 
Otherwise, if gas is limited, the system does not have much 
fuel to grow its BH/bulge mass, and its position on the Mpgg- 
Mestar diagram remains roughly the same. In this case, since 
the z = 0 position is near the Mpy-Mstar relation, the cur- 
rent position should also be near the relation (Fig. | right). 

In this work, we assess the gas content of bulge-dominated 
galaxies using ALMA observations. Our ALMA observa- 
tion targeted CO(2-1) or CO(3-2) for 7 bulge-dominated 
galaxies. For our targets, the low-J transitions are the best- 
possible lines accessible by ALMA to estimate Mgas, as they 
are not subject to strong assumptions about the CO spectral 
line energy distribution (SLED). Based on the ALMA data, 
we constrain the molecular-gas content and discuss the im- 
plications for BH-bulge coevolution. We note that each of 
the 7 individual systems should follow the BHAR-SFR rela- 
tion (see Footnote 1), and thus they are suitable to test the 
two scenarios laid out in Fig. 1. 

The structure of this paper is as follows. In $2, we reduce 
and analyze the ALMA data. We also compile existing mul- 
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tiwavelength data and perform a spectral energy distribution 
(SED) analysis. In $3, we discuss the physical implications 
of our results. We summarize our results and discuss future 
prospects in $4. 

Throughout this paper, we assume a cosmology with Ho = 
70 km s^! Mpc-!, Qj = 0.3, and Q4 = 0.7. We adopt a 
Chabrier initial mass function (IMF; Chabrier 2003). Quoted 
uncertainties are at the 1a (68%) confidence level. The word 
"gas" in this work specifically means cold molecular gas 
(mostly Hə and He), unless otherwise stated. 


2. DATA AND ANALYSIS 
2.1. Targets and observations 


Our ALMA Cycle 7 program (2019.1.00678.S, PI: G. 
Yang) targeted seven bulge-dominated star-forming galax- 
ies in the GOODS-South field, which has excellent multi- 
wavelength coverage (see $2.5). The targets were selected 
from the bulge-dominated sample in Yang et al. (2019), 
and their basic properties are summarized in Table. 1. The 
bulge-dominated sample is classified using machine-learning 
morphological measurements (Huertas-Company et al. 2015) 
based on HST H-band imaging (Grogin et al. 2011; Koeke- 
moer et al. 2011). The machine-learning approach is efficient 
and reliable and has been successfully applied to other fields 
beyond GOODS-South (e.g., Huertas-Company et al. 2015; 
Ni et al. 2021). This targeting approach is designed to select 
pure spheroidal galaxies that do not have a disk component 
(see, e.g., Fig. 2 of Yang et al. 2019), and we discuss conse- 
quences of missing a disk component in $3.1. It is essential 
to focus on bulge-dominated galaxies, as the BHAR-SFR re- 
lation does not exist among other (e.g., disky and irregular) 
galaxies (Yang et al. 2019). Consistently, in the local uni- 
verse, BH masses are only tightly correlated with the masses 
of bulges rather than disks (e.g., Kormendy & Ho 2013). 

Our ALMA targets are selected to have secure optical 
spectroscopic-redshift (spec-z) measurements from the liter- 
ature, and the redshifts span z = 0.5-2.5 (see Table. 1). The 
spec-z measurements are necessary to locate the observed- 
frame CO frequencies and target them with ALMA. The sam- 
ple is also selected to have Spitzer/MIPS and/or Herschel 
(S/N > 3) detections (82.5), so that they are likely at the 
early star-forming phase after bulge formation ($1). The Her- 
schel flux uncertainties are based on Monte Carlo simulations 
that account for both instrumental and confusion noise (El- 
baz et al. 2011). We are primarily interested in star-forming 
rather than quiescent bulges. This is because, for quiescent 
systems, both bulge and BH growth have essentially ceased 
according to the BHAR-bulge SFR relation (Yang et al. 2019) 
and the system is already in place on the Mpy-Mstar re- 


lation.? Five (out of seven) targets are classified as X-ray 


AGNs in Luo et al. (2017). This prevalence of AGNs is 
naturally expected from the BHAR-bulge SFR relation, as 
our targets are selected to be star-forming bulge-dominated 
galaxies. 

The ALMA data were taken on January 14 and 24, 2020, 
with the on-site perceptible water vapor (PWV) ranging from 
4.6 mm to 6.6 mm. We target CO(2-1) and CO(3-2) for 
z « 2 and z > 2 sources in our sample, respectively. These 
lines are within the frequency ranges covered ALMA bands 
3 and 4. Among the four spectral windows (spw), one spw 
was placed centered on the line, and the other three covered 
continuum frequencies. These observations were performed 
on a 12-m array configuration of C-3 (maximum baseline — 
0.5 km), resulting in angular resolutions of 1.1—2.3 arcsec. 
The on-source exposure time was 24—28 minutes, reaching a 


1c continuum sensitivity of ~ 0.02—0.03 mJy beam !. 


2.2. Data reduction 


We requested and downloaded the calibrated Measurement 
Sets (MS) for our observations using the online Science 
Ready Data Products (SRDP) service. These MS data are 
produced by the ALMA Pipeline v6.1.2-7. We use the Com- 
mon Astronomy Software Applications (CASA) v6.2 package 
to further reduce these MS data. 

We employ the “tclean” function in CASA to produce the 
primary-beam (PB) corrected continuum images from the 
visibility data. We include all available spw (masking poten- 
tial line channels within +500 km s~! around the line cen- 
ter). We set the output image size to 20" x 20" with a pixel 
scale of 0.2". We then apply the “imfit” function (CASA) 
to these images to determine the source continuum position. 
Five sources have converged fits but none of them is signifi- 
cant (S/N > 3), indicating the continuum emission is weak. 
We present a quantitative analysis of the continuum fluxes in 
824. 

For the line analysis, we first utilize tclean to transform 
the visibility data to an image cube for each source. We 
use the line spw and set the output spectral channel width to 
50 km s~! (dataset native resolution ~ 1-3 km s71), while 
the spaxel size and pixel scale are the same as above. Using 
“imcollapse” (CASA), we then collapse the cube along the 
frequency axis including the line frequency (inferred from 
the optical redshift) +300 km s^!. Based on the resulting 
line image, we perform source detection with “imfit” (CASA), 


? The program actually observed 10 targets in total. These 10 targets were 
originally selected using an earlier IR catalog, PEP (Lutz et al. 2011). How- 
ever, the most recent catalog (Barro et al. 2019; $2.5), which carefully ad- 
dresses source confusion, indicates 3 targets are actually quiescent, being 
undetected in the IR. Therefore, we only focus on the rest of the targets, 
i.e., 7 IR star-forming galaxies. 


3 https://data.nrao.edu/portal/#/ 
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Scenario 1: the BHAR-SFR relation 
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Scenario 2: the BHAR-SFR relation 
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Figure 1. Schematic diagram of the two scenarios tested in this work. The blue and red stars represent the bulge-formation redshift and 
z = 0, respectively. The black arrows indicate the evolution paths. The dashed line represents the observed BH-bulge mass relation in the local 
universe. The determining factor between the two scenarios is the amount of cold gas available after bulge formation (see §1 for details). 


Table 1. Basic source properties 


ID RA DEC Z z ref. log Mstar log SFR log Lx Tgrow tua 

(1) (2) (3) (4) (5) (6) (7) (8) (9) (10) 
528 53.113468  —27.933294 1.089 Cooper etal. (2012) 11.28 +0.10 1.58+0.04 41.754+0.16 509+1.25 5.42 
6278 53.060116 —27.852997 1.540 Suh et al. (2015) 10.97+0.05 1.42+0.23 43.00+0.07 3.524 1.88 4.11 
23845 53.097649  —27.715282 2.142 Coil et al. (2015) 10.96 +0.02 2.05+0.02 43.412- 0.09 0.80+0.06 3.02 
24210 53.071419 —27.717581 0.566 | Cooperetal.(2012) 10.58 +0.05 1.89+0.02 42.87 0.05 . 0.49 0.06 7.98 
24682 53.104096 —27.683758 0.732 Cooper etal. (2012) 10.78 +0.02 0.58+0.04 42.98 +0.05 15.73 1.62 6.99 
25573 53.139187  —27.694145 1.044 Vanzellaet al. (2008) 11.00 +0.02 1.14 Æ 0.04 — 7.234 0.76 5.58 
25998 53.137573 —27.700104 2.453 Barro et al. (2013) 11.05 +0.07 2.33+0.06 43.68-E 0.08 0.53 40.11 2.63 


NOTE. — (1) Identification in the CANDELS catalog (Guo et al. 2013). (2) & (3) CANDELS J2000 coordinates. (4) Optical spectroscopic 
redshift. (5) Redshift reference. (6) & (7) Logarithmic stellar mass (Mo) and star formation rate (Mo yr~') from our SED modelling (see 


82.5). (8) Intrinsic X-ray luminosity (erg s71) based on the absorption-corrected X-ray flux (see $2.5). 


€ o» 


means X-ray undetected. We note 


that Lx only reflects instantaneous AGN activity, not long-term average BH growth (see $1). (9) The bulge stellar growth timescale (Gyr) as 


defined in Eq. 5. (10) The Universe's age (Gyr) at the source's redshift. 


which searches for a Gaussian-shaped source. If a source is 
detected with signal-to-noise (S/N) > 3, we adopt the imfit 
source position for spectrum-extraction below; otherwise, we 
adopt the CANDELS position (Table 1). 

The empirical choice of collapsing the data cube within 
+300 km s~! is to cover most of the line signal, as CO line 
width is typically < 600 km s^! in the literature (e.g., Fre- 
undlich et al. 2019; Shangguan et al. 2020). Choosing an 
even larger width could dilute the S/N. We also test different 
ranges other than +300 km s~! (see $2.3) and find similar 
results. Admittedly, if a line center has a large shift greater 
than 300 km s-!, our method could miss it. However, from 


our extracted spectra, we do not find such strong shifts (see 
$2.3). 


2.3. Line-flux measurements 


We extract the CO spectrum for each source with CASA's 
"specflux", utilizing a circular aperture with an area of 
2xbeam area. Here, we choose a fixed aperture instead of 
an adaptive aperture based on, e.g., contours. This is be- 
cause our targets are generally faint and affected by random 
noise significantly. The latter adaptive approach, which fa- 
vors "positive" noise and avoids "negative", may potentially 
introduce a positive bias. This bias could lead to S/N over- 
estimation and even false detections. The resulting spectra 
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CANDELS 528, Ico = 706 + 88 mJy km s-! 
1.0} : 


f, [my] 


f, [mJy] 


f, [mJy] 


f, [mJy] 


-1.0 


f, [mJy] 


-0.5 


CANDELS 25998, Ico = 602+ 95 mJy km s- 
—1000 —750 —500 —250 0 250 500 750 1000 


velocity [km s71] 


Figure 2. Left: HST H-band 7" x 7” cutouts with contours of CO emission. The contours are at the 2, 3, 5, and 8 sigma levels. The beam 
profile is displayed at the bottom-left corner. There is no contour for CANDELS 23845 due to the weak signal of its CO-line map. Right: CO 
spectrum for each source. The red dashed curve represents the best-fit Gaussian model (only displayed for CO detected sources). The horizontal 
and vertical dotted lines indicate zero flux and velocity, respectively. The orange shaded region indicates the integrated velocity range for the 
line-flux measurement. The measured line flux and its uncertainty (or upper limit) is labeled. We consider the line as detected if S/N 3; we 
adopt the 3c uncertainties as upper limits for undetected sources (see Table 2). 
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are displayed in Fig. 2. From these spectra, there do not ap- 
pear to be any significant line signals beyond 300 km s71, 
supporting our source-detection method in $2.2. 

We then fit the CO spectrum for each source with CASA's 
"specfit" employing a Gaussian model. Here, we do not in- 
clude a continuum component, because the continua are in- 
significant for all of our sources (see $2.4). To avoid false 
detections, we require that the Gaussian center (Ucent) is 
within +200 km s^!, where the velocity zeropoint corre- 
sponds to the line frequency inferred from the optical red- 
shift (and hereafter). The fits converge for 6 (out of 7) 
sources. We obtain their best-fit Vcent and full width half 
maximum (FWHM) from the specfit result. For each of the 6 
sources, we estimate the velocity-integrated line flux (units: 
Jy km s^!) by integrating the CO spectra over the range of 
(Ucent -FWHM, Veent +KWHM). For the other source (CAN- 
DELS 23845), we perform the integral over (—300 km s^ !, 
+300 km s71). 

To estimate the line-flux uncertainty for each source, we 
randomly place ~ 20 non-overlapping apertures around each 
source (avoiding its spectrum-extraction region). For each 
aperture, we then extract the spectrum and measure the 
velocity-integrated flux in the same way as above. We calcu- 
late the standard deviation of these fluxes. We repeat this pro- 
cess 100 times obtaining 100 standard-deviation values, and 
then calculate the median of these values. Finally, we adopt 
the median as the CO flux lo uncertainty for each source. 
If a source has CO-flux S/N» 3, we consider the line as de- 
tected, and 4 out of 7 sources have CO detected. For the other 
3 sources, we adopt 3x CO-flux errors as their upper limits. 
In $2.2, we collapse the data cubes over +300 km s-! for an 
initial source search. We also test collapsing at other velocity 
ranges from +100 km s^! to +500 km s^! but do not find 
additional CO detections beyond the 4 sources. Therefore, 
we conclude that our results are not sensitive to the collapsed 
velocity range. 

Since the line flux above is measured within an aperture 
(“aperture flux”, hereafter), we need to do an aperture cor- 
rection to account for the emission outside the aperture. To 
perform this task, we use a large circular aperture with an 
area of 4x beam area to estimate the “total” line flux for our 
highest S/N source (CANDELS 528 with S/N ~ 8). We then 
divide this total flux by the aperture flux and adopt the result 
as our aperture-correction factor (— 1.12). We multiply the 
line flux/error (or upper limit) by this aperture-correction fac- 
tor. We also test estimating the correction factor based on the 
other 3 detected sources. The resulting corrections are 1.10— 
1.28, similar to our adopted value (1.12). Our adopted cor- 
rection is based on the highest S/N source, and thus should 
be the most reliable. We note that our main results ($3) are 
not sensitive to the choice of the aperture-correction factor, 
as all the values (from 1.10 to 1.28) are relatively small. Ta- 


ble 2 lists the final CO fluxes, errors, and upper limits for our 
sources. For the four CO-detected sources, we present spa- 
tial analyses of their CO emission in Appendix A. In brief, 
the angular resolutions of our ALMA observations are not 
sufficient to well resolve the detected CO emission. 

From the measurements of CO line fluxes above, we es- 
timate CO luminosities following Solomon & Vanden Bout 
(2005), i.e., 

LoolJ, J —1) = 3.25 x 10 Icora D2(1--z) 3, (1) 


obs 


where /co is the velocity-integrated line velocity in 
Jy km s^ ry is the line frequency in the ob- 
served frame; Dy is the luminosity distance in Mpc; 
Loo(J, J — 1) is in units of K km s^! pc?. Assum- 
ing To1 = Leo (2, 1)/Loo (1, 0) = 0.8 and T31 = 
Liao (3, 2)/Leo (1,0) = 0.5 (e.g., Saintonge et al. 2017; 
Lamperti et al. 2020), we can convert the observed 
Loo (J, J — 1) (J = 3 or 2) to Log (1,0) (Leo hereafter). 
Finally, we estimate the molecular mass from 


Mas E acoLco: (2) 


where we assume the conversion factor oco = 
4.3 Mc (K km s^! pc?)~! (we use these units for aco, 
hereafter), a typical value adopted in the literature (e.g., 
Bolatto et al. 2013; Carilli & Walter 2013). This adopted 
conversion factor includes the contribution from Helium and 
metals. We discuss the effects of this aco assumption in 
§3.1. Table 2 lists the resulting Loo and Meas. 


2.4. Continuum emission 


In §2.2, we performed a source search on the ALMA con- 
tinuum images but did not find any significant detections. 
Therefore, we adopt the CO line position (if available) or the 
CANDELS position for the continuum measurements. As 
for the line-flux extraction, we also employ a circular aper- 
ture with area = 2xbeam area to measure the continuum 
flux. We estimate the noise by randomly placing apertures 
and calculating the standard deviation of the resulting con- 
tinuum fluxes (a similar procedure as in §2.3). None of our 
sources has S/N> 3 continuum fluxes, consistent with the 
result in $2.2. As in $2.3, we also adopt 3 x noise as the con- 
tinuum upper limit for each source (listed in Table 2). We 
then convert these flux upper limits to Ly 850 um (rest-frame 
850 um luminosity) upper limits, assuming the Rayleigh- 
Jeans law (L, c v?) for K corrections. Finally, we con- 
strain the continuum-based gas masses using the scaling rela- 
tion, Moon* = Ly 850 ,m/ (1.01 x 10% erg s! Hz-! M5"), 


gas 


from Scoville et al. (2016, 2017).4 


^ The original scaling factor was 6.7 x 101? erg s71 Hz-! Ms! assuming 
aco = 6.5 (Scoville et al. 2016). Here, we modify the scaling factor so 
that the relation becomes consistent with our assumed aco = 4.3 (82.3). 
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Table 2. ALMA results 


ID J Amp. Ucent Av Ico log Megas Tdep u SS log M. oo 

(D Q (3) (4) (5) (6) (7) (8) (9) (10) (11) 
528 2 1.814 0.22 35 4224 4044056 706-88 10.78 +0.05 158+0.24 0.31+0.08 < 79 « 10.75 
6278 2 - - -— « 155 « 10.41 < 0.96 «0.27 <70 « 10.98 
23845 3 - « 212 « 10.65 « 0.40 «0.49 «80 « 10.95 
24210 2 1.39+0.24 —42+ 26 303 +62 409+ 82 9.97+0.09 0.12+0.02 0.24+0.06 <93 < 10.25 
24682 2 0.9240.24 —62 +28 223467 193+ 63 9.87 +0.14 1.9340.65 0.12+0.04 <95 < 10.49 
25573 2 - - - « 151 « 10.07 « 0.86 «0.42 <77 « 10.70 
25998 3 0.92+0.15 —33+51 618+119 602495 11.21 +0.07 0.76+0.16 143+0.32 < 88 < 11.10 


NOTE. — (1) Identification in the CANDELS catalog. (2) Targeted CO transition (J — J — 1). (3), (4), & (5) Gaussian amplitude (mJy), 


central velocity (km s7 +; relative to that from optical redshift), and FWHM (km s+) from the fit of the CO spectrum. “—” indicates 
CO-undetected (S/N < 3). (6) Velocity-integrated CO flux or 3c upper limit (if S/N < 3) in mJy km s^ !. (7), (8) & (9) Logarithmic gas mass 
(Mo) or 3c upper limit (inferred from Lgo; see $2.3), gas-depletion timescale (Gyr), and gas-to-stellar mass ratio. (10) & (11) 3c upper limit 


of continuum flux (jJy) and corresponding logarithmic gas mass (Mo). 


The resulting MID upper limits are listed in Table 2. 
For 3 out of the 4 CO-detected sources (except CAN- 
DELS 25998), the CO-based masses are consistent with 
the continuum-based constraints. For CANDELS 25998, 
the CO-based mass is slightly higher than the continuum- 
based upper limit by ~ 0.1 dex (see $3.1 for more dis- 
cussion of this object). For the 3 CO-undetected sources, 
the CO-based mass constraints are tighter than continuum- 
based ones. Therefore, our scientific discussions are primar- 
ily based on the CO measurements in $3, unless otherwise 
stated. 


2.5. Multiwavelength Spectral Energy Distributions 


We employ CIGALE v2022.0 (Roehlly et al. 2014; Boquien 
et al. 2019; Yang et al. 2020, 2022) to perform SED mod- 
elling for our ALMA sources. We compile multiwavelength 
broad-band photometric data from the Rainbow Cosmolog- 
ical Surveys Database. These data include 23 bands from 
U to SPIRE 500um (Guo et al. 2013; Barro et al. 2019). 
We also compile X-ray photometry from the 7 Ms Chan- 
dra Deep Field-South (CDF-S) catalog (Luo et al. 2017). 
Six sources (all except CANDELS 25573) have X-ray de- 
tections. Among them, five are detected in the hard band 
(2-7 keV) and we adopt their hard-band fluxes, which are 
less affected by obscuration compared to the full-band (0.5— 
7 keV) and soft-band (0.5-2 keV) fluxes. Since CIGALE 
requires obscuration-corrected X-ray fluxes as input (Yang 
et al. 2020, 2022), we apply obscuration corrections to these 
adopted fluxes. The corrections are estimated based on the 
absorption column densities (Np) from Luo et al. (2017) and 


5 http://arcoirix.cab.inta-csic.es//Rainbow. navigator. public/ 


PIMMS.Ó One source (CANDELS 528) is only detected in the 
soft band, and we adopt the soft-band flux." 

The model parameters in CIGALE are listed in Table 3. 
We adopt a delayed-r model (sfhdelayed in CIGALE) 
for star formation history (SFH) and a Bruzual & Charlot 
(2003) model (bc03 in CIGALE) for a simple stellar popula- 
tion (SSP). We note that, when fitting the data, CIGALE auto- 
matically excludes the unphysical models in which the stellar 
age is older than the Universe's age. In bc03, we assume a 
Chabrier (2003) IMF and a solar metallicity (Z — 0.02). We 
also include nebular emission (nebular in CIGALE) with 
default settings. We adopt the Charlot & Fall (2000) model 
(dustatt_modified_CF00 in CIGALE) for stellar atten- 
uation. We allow the V -band attenuation varying from 0.2 to 
3 mag, while leaving other parameters (uv, nism, and npc) 
at the default values. We use the Dale et al. (2014) model 
(dale2014 in CIGALE) for galactic dust emission. We al- 
low three values for radiation-field slope, i.e., 1.5, 2, and 2.5. 
We use the SKIRTOR model (Stalevski et al. 2012, 2016, 
Skirtor2016 in CIGALE) for the AGN UV-to-IR emis- 
sion. We allow fracagn (fractional AGN IR luminosity) to 
vary between 0 and 0.99. We set ÜAcw (viewing angle) to 
30? and 70°, which are typical values for type 1 and type 2 
AGNS, respectively (e.g., Yang et al. 2020; Ramos Padilla 
et al. 2021). Finally, we include the xray module in CIGALE 
to account for AGN/galaxy X-ray emission, and leave the re- 
lated parameters at the default value(s). The CIGALE configu- 
rations above lead to a total of 27,243,160 models (3,891,880 


6 https://cxc.harvard.edu/toolkit/pimms.jsp 


7 Obscuration correction is not possible for this source due to the single-band 
detection. However, we do not expect the obscuration is strong, because 
otherwise the detected band would likely be the hard band rather than the 
soft band. 
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per source). In addition to the photometry above, we also in- 
clude the ALMA continuum upper limits ($2.4) in the fits. 

We run CIGALE with the configurations above. For 4 
out of 7 sources, the fit quality is acceptable with reduced 
x? < 2 (see Fig. 3 for their best-fit SEDs). We adopt the 
Bayesian output of Mstar and SFR (Table. 1)3 However, 
for CANDELS 23845, 24210, and 25998, the observed far- 
IR fluxes are systematically higher than the model values, 
leading to large reduced x? zz 4—5. This “far-IR excess” 
might be caused by a dust-enshrouded stellar population that 
is strongly attenuated at shorter wavelengths (e.g., Hodge 
et al. 2016; Buat et al. 2019). We note that CANDELS 23845, 
24210, and 25998 are actually similar to the sample of Buat 
et al. (2019), i.e., they are all massive (log(Mstar/Mo 2 
10.5) dust-rich (log(Lin/Lco) Z 12) galaxies. The Buat 
et al. (2019) sample's continuum emission is detected by 
ALMA, but ours is not ($2.4). We consider this difference 
is due to the different ALMA bands being used (band 6 vs. 
band 3) and that dust emission is stronger at shorter ALMA 
wavelengths due to the Rayleigh-Jeans law (e.g., Fig. 2 of 
Scoville et al. 2016). 

For these three sources, we include an additional ad hoc 
modified (flexible emissivity) black body model (mbb in 
CIGALE) to account for this far-IR excess (Table 3) and re- 
run CIGALE. This mbb component is designed to model the 
far-IR emission from the dust-enshrouded stellar population 
as discussed above. To account for possible temperature vari- 
ance of this hidden population, we allow a cold model (50 K) 
and a hot model (100 K) in the mbb (see Table 3). We al- 
low the luminosity ratio, Dmpp/Laale2014, to vary from 0 to 
50. Indeed, the new fits have improved fit quality (reduced 
X? = 1-2; see Fig. 3) compared to the previous fits. After 
adding the new component, the Akaike information criterion 
(AIC; Akaike 1974) is reduced by > 10, indicating the im- 
provements with mbb are statistically significant (e.g., Burn- 
ham & Anderson 2002). 

For the new fits, we do not adopt the SFRs directly from the 
CIGALE output, because the CIGALE SFRs do not account for 
the mbb contribution which dominates the galaxy IR lumi- 
nosity (Li,gai)). We estimate SFRs based on the fitted total 
galaxy IR luminosity, excluding the AGN contribution (Ken- 
nicutt 1998; Salim et al. 2007), i.e., 


SFR = 1.09 x 1071 Ei gal (3) 


where SFR and Ln, gal are both in solar units. We still adopt 
the output Mstar, assuming that the stellar mass is dominated 
by the main stellar population rather than the hidden popu- 
lation. This assumption is justifiable, because the NIR data 


8 The Bayesian output is a probability-weighted value that considers all mod- 
els available. It is thus more robust than the best-fit output, which is based 
on a single model with minimum x?. 


(rest-frame ~ lum), a robust Mstar indicator, can be fitted 
well with the main stellar population alone (see Fig. 3). We 
note that our main conclusions are unaffected, even if we 
miss some Mstar contributed by the dust-enshrouded stellar 
population (see $3.1). 


We note that another approach to account for the far-IR 


excess is to adopt a flat dust-attenuation curve, which leads 
to significant attenuation at near-IR wavelengths (e.g., Buat 
et al. 2019). This method effectively assumes that the far- 
IR excess comes from the strong attenuation in the near-IR. 
We also test this approach by allowing shallower ISM atten- 
uation slopes (nism) ranging from —0.3 to —0.7 (disabling 
mbb). Indeed, the fit quality has been improved (reduced 
x? = 2-2.5) compared to the original fits (reduced x? ~ 4— 
5). The resulting SFRs (also estimated from Lrg, ga1) are sim- 
ilar to those from the mbb fits, with differences < 0.1 dex, 
but the Mstar are systematically higher by ~ 0.3-0.6 dex.? 
These effects are similar to what Buat et al. (2019) found. In 
$3.1, we discuss the effects to our main conclusions, if the 
flat-curve approach is adopted instead of the mbb one. 


The Mstar and SFRs for all sources are listed in Table 1. 


Following Ni et al. (2021), we classify a source as a star- 
forming galaxy if it satisfies SFR > SFRys/10! 4, where 
SFRys is the main-sequence SFR as a function of Mstar 
and redshift from Whitaker et al. (2012). All of our 7 
sources are classified as SF galaxies. From Table 1, CAN- 
DELS 6278 has the largest SFR uncertainty (0.23 dex com- 
pared to < 0.05 dex for the other sources). This is un- 
derstandable, as CANDELS 6278 is the only object without 
Herschel detections in our sample (see Fig. 3). This level 
of SFR uncertainty (0.23 dex) without Herschel is realistic 
based on previous studies with CIGALE (e.g., Mountrichas 
et al. 2021). 


3. DISCUSSION 
3.1. Implications for BH-bulge coevolution 


With the measurements of galaxy properties in $2, we now 


discuss the implications for BH-bulge coevolution. First, we 
define a gas-depletion timescale as 


Meas 
ep ==) 4 
Tdep SFR ( ) 
and bulge stellar growth timescale as 
— Mstar (5) 
"erc ep 


? Based on these alternative estimations of SFRs and Mstar, the sources are 
still classified as star-forming galaxies. 

10 The error is calculated by CIGALE as the standard deviation of the marginal- 
ized SFR probability distribution, accounting for all available physical 
models (Boquien et al. 2019). 
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Table 3. CIGALE model parameters 


Module Parameter Symbol Values 
Star formation history Stellar e-folding time Tatar 0.5, 1, 2, 3, 4, 5 Gyr 
sfhdelayed Stellar age tstar 0.2, 0.5, 1, 2, 3, 4, 5 Gyr 
Simple stellar population Initial mass function — Chabrier (2003) 
bc03 Metallicity Z 0.02 

V -band attenuation in the ISM ADM 0.2-1 (step 0.1), 1.5, 2, 2.5, 3 mag 
Dust attenuation ADM CAD + ABE) uv 0.44 
dustatt.modified.CF00 Slope of the ISM attenuation NISM —0.7 

Slope of the birth-cloud attenuation nBC —1.3 

: Pa Slope in dMaus x U7 “test dU Odust 1.5, 2, 2.5 

oo Ratio of Lipp and Laaie2014 Ed 0, 1,2, 5, 10, 20, 50 
mbb? Temperature of mbb Tabb 50, 100 K 

Emissivity of mbb B 1.5 

AGN contribution to IR luminosity fracaan 0—0.9 (step 0.1), 0.99 


AGN (UV-to-IR) emission ae 
skirtor2016 Viewing angle 


Polar-dust color excess 


ÜAGN 30°, 70° 
E(B—V)pp 0.03, 0.2, 0.4 


X-ray emission 


xray AGN X-ray angle coefficients 


Maximum deviation from the &ox-L, 25004A relation | Adox| max 0.2 


(a1, a2) (0.5, 0) 


NOTE. — For parameters not listed here, we use the default values. (a) The mbb module is only used for CANDELS 23845, 24210, and 


25998 to improve their fit quality (see §2.5). 


From this definition, Tgrow represents the timescale needed to 
double the stellar mass in the future, given the current SFR. 
Denoting the galaxy’s redshift corresponds to a Universe’s 
age of to and assuming the SFR keeps constant with suffi- 
cient gas supply, then at cosmic time of t < to + Tgrow, 
the stellar mass is predominantly assembled before to; at 
t > to Tgrow, the mass is predominantly assembled after to. 
This stellar growth timescale is also a proxy for BH growth 
timescale due to 


Mstar Mstar /300 "m Mpggn 


Teow = SFR ^ SFR/300  BHAR 


= Tgrow,BH; (6) 


where we apply the long-term average BHAR-SFR relation 
(BHARz SFR/300; see $1) and assume the local BH-bulge 
mass relation (Mpy 7 Metar/300). We note that the system 
might not follow the local BH-bulge mass relation, and thus 
Eq. 6 is only an approximation. If Mpun is higher (lower) 
than the local BH-bulge mass relation, then Tgrow,BH should 
also be longer (shorter) than Tgrow. The values of 75,4, and 
Tdep for our targets are listed in Tables 1 and 2, respectively. 

By comparing the two timescales of Tacp and Tgrow, we 
can gain insight the BH-bulge coevolution. If 7465 > Tgrow, 
then the gas content can last sufficiently long for signifi- 
cant BH/bulge growth, as the BH/bulge mass predominantly 
forms during the bulge phase where the BHAR-SFR relation 
applies (Scenario 1 in Fig. 1). Otherwise, the gas is depleted 
quickly before significant BH/bulge growth and the system 
cannot significantly change its position on the Mpy-Mstar 


diagram (Scenario 2 in Fig. 1). Fig. 4 compares Tqep versus 
Tgrow for our ALMA targets. Six (out of seven) SF sources 
have Tgep < Tgrow, With the only exception being CAN- 
DELS 25998. For these 6 sources, their Tgep is more than 
2 times shorter than Tgrow.- 

Our estimates of Mgas assume a typical CO-to-gas con- 
version factor of aco = 4.3. Some studies suggest a lower 
value of aco = 0.8 for compact galaxies (e.g., Bolatto et al. 
2013; Carilli & Walter 2013). Considering that our bulge- 
dominated galaxies are generally compact (e.g., Ni et al. 
2019, 2021), we also plot the result for the assumption of 
aco = 0.8 in Fig. 4. Under this alternative assumption, 
all seven sources, including CANDELS 25998, have Tgep at 
least 4 times shorter than Tgrow- 

Tdep and Tgrow are measures of what can happen in the fu- 
ture life of the galaxy assuming the current SFR. But SFR 
is unlikely to be constant over cosmic time. From hydrody- 
namical simulations, SFR can be highly variable on relatively 
short timescales of ~ 10 Myr (e.g., Flores Velázquez et al. 
2021). However, from Eqs. 4 and 5, both 7445 and Tgrow are 
inversely proportional to SFR, and thus their ratio is indepen- 
dent of SFR, i.e., 

Tdep _ Meas - 

Terow = Mara =H, (7) 
and u is also known as the gas fraction of the galaxy. To 
compare Tgep With Tgrow, it is equivalent to compare u with 
unity. The advantage of ju is that it is not affected by the SFR 
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Figure 3. SED fits of our ALMA sources using CIGALE. The black curve represents the best-fit SED model, while the blue, orange, and red 
curves indicate AGN, galaxy, and mbb (if present) components, respectively. The redshift and reduced X? are labeled on each panel. 
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variability mentioned above. Aside from this mathematical 
feature, y = Meas /Mstar also has a useful physical mean- 
ing: it represents the amount of mass the BH/bulge can po- 
tentially grow in the future compared to the current mass. If 
u is above unity, the "potential" mass will be dominant, sup- 
porting Scenario 1 in Fig. 1; otherwise, the "current" mass 
will be dominant, supporting Scenario 2 in Fig. 1. 

Fig. 5 shows u versus redshift for our ALMA targets, un- 
der the assumption of both aco = 4.3 and aco = 0.8. As 
expected, all sources have u < 1 except CANDELS 25998 
(for aco = 4.3). Considering the upper-limit data points, 
the median values of js are only < 0.26 (for aco = 4.3) and 
< 0.048 (for aco = 0.8), significantly smaller than unity 
under both aco assumptions. Another source of systematic 
uncertainties is the assumed CO line ratios. Our assumed 
values (rə1 = 0.8 and r3; = 0.5) are typical among star- 
forming galaxies (e.g., Saintonge et al. 2017; Lamperti et al. 
2020). Some observations suggest that AGN hosts (the case 
for many of our targets; Fig. 3) tend to have rg; and r31 
higher than normal galaxies (e.g., Kirkpatrick et al. 2019). 
We note that adopting higher ro; and r3; would decrease the 
final gas masses and u, and thereby strengthen our main re- 
sult above. 

In the discussion above, we assume that total Mstar œ% 
bulge Mstar, because our targets are bulge-dominated from 
HST imaging ($2.1). However, it is still possible that a minor 
disk component is missed, which could account for « 40% 
of the total mass (e.g., Huertas-Company et al. 2016). In 
the most extreme case (40% mass from a missed disk com- 
ponent), our bulge Mstar is overestimated by a factor of 
1/(1—0.4) = 1.67, and thus pz is underestimated by the same 
factor according to Eq. 7. Even this level of  underestima- 
tion is not strong enough to overturn our main result, as most 
sources (except CANDELS 25998 assuming aco = 4.3) 
will still have u < 1 after correcting for the underestimation 
(see Table. 2). 

In Fig. 5, we include two star-forming bulge-dominated 
sources also in the GOODS-South field from Barro et al. 
(2016).! One source is CANDELS 21662 at z = 2.18. 
The other is CANDELS 25998 at z = 2.45, which is also in 
our sample. The Barro et al. (2016) Mgas values were from 
Maust assuming a gas-to-dust ratio of 100, and their Maust 
estimation is based on SED modeling of ALMA continuum 
and other IR data. For CANDELS 25998, the Barro et al. 
(2016) Megas (1019? Mo) is consistent with our continuum- 
based Mg, constraint («101 -! Mo), although the two 
measurements are based on different frequencies and as- 
sumptions. The two u values from Barro et al. (2016) are 


!! Barro et al. (2016) have six star-forming galaxies in total, but only two of 
them are bulge-dominated according to the morphological classifications 


of Huertas-Company et al. (2015). 


both below unity, consistent with our sample. For CAN- 
DELS 25998, the Barro et al. (2016) Mstar estimation is sim- 
ilar to ours (1011-97 vs. 1011-05 Mọ), and thus the cause of 
the u difference between Barro et al. (2016) and ours (the 
rightmost points in Fig. 5) is mostly related to the Mgas mea- 
surements. The Barro et al. (2016) u is between our values 
based on aco = 0.8 and aco = 4.3, suggesting that the 
intrinsic @co for this source is within 0.8—4.3. But we note 
that the uncertainties of Barro et al. (2016) yu are relatively 
large, and both of our u values are consistent with their mea- 
surements at a 2c confidence level. 

Our results above suggest that, without gas replenishment 
(see $3.2), the cold molecular gas of the SF bulge-dominated 
galaxies will be depleted before significant BH/bulge growth. 
Their Mpu/Mstar ratios remain largely unchanged during 
the bulge-evolution phase until z = 0 (see $1). Therefore, 
it is likely that the Mpy-Mstar relation has already formed 
at the beginning of the bulge phase and they maintain this 
relation until z = 0 (i.e., Scenario 2 in Fig. 1). The de- 
tailed formation mechanisms are unknown, subject to inves- 
tigations from both theoretical and observational approaches 
(e.g., Huertas-Company et al. 2018; Hopkins et al. 2022). 
The BHAR-bulge SFR relation has the role of maintaining 
the BH-bulge mass correlation, but not creating it. We cau- 
tion that our sample is limited to z = 0.5-2.5. It is possible 
that this conclusion might change at higher redshift, consid- 
ering that cold gas tends to be more abundant toward the early 
universe. Also, from Fig. 5, the z > 2 sources tend to have 
higher u than sources at lower redshifts, although our sample 
size not sufficiently large to robustly test this trend. 

In our SED fits ($2.5), we employ an ad hoc mbb compo- 
nent to account for the FIR excess in three sources. We might 
underestimate Mstar as we do not account for the mass con- 
tribution from the hidden stellar population. Therefore, the 
actual p values might be even lower than our case, further 
strengthen our main conclusion. The FIR excess can also be 
addressed by adopting a flat attenuation curve. We note that 
this approach would also lead to higher Mstar estimations 
(see $2.5) and thereby lower u. In summary, we consider 
that our main conclusion is not affected qualitatively by the 
details of the SED-fitting procedure. 


32. Gas replenishment 


The discussion in $3.1 assumes that the gas content in our 
bulge-dominated galaxies is primarily consumed by star for- 
mation without replenishment. However, if gas replenish- 
ment is common and supplies Mgas by a typical factor of 
= 4 (as median u < 0.26 for aco = 4.3 in 83.1), then the 
BH/bulge growth can sustain for a longer timescale than our 
estimated Tgep. 

This possibility can be qualitatively investigated by study- 
ing the fraction of star-forming galaxies, as common gas re- 
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Figure 4. 74.5 versus Tgrow. The red stars and squares are values 
estimated using aco = 4.3 (default) and aco = 0.8, respectively. 
The downward-pointing arrows indicate 3c upper limits. The black 
dashed lines indicate Taep = Tgrow. Except for CANDELS 25998 
(aco = 4.3), other data points are all below the Tacp = Tgrow line. 
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Figure 5. u = Mgas/Mstar versus redshift. The stars and squares 
are values estimated using aco = 4.3 (default) and aco = 0.8, re- 
spectively. The downward-pointing arrows indicate 3c upper limits. 
The dashed line indicates u = 1. The two blue data points repre- 
sent two bulge-dominated galaxies at z — 2—2.5 from Barro et al. 
(2016). The gas content of the bulge-dominated galaxies is gener- 
ally low (u < 1 in most cases), insufficient to support significant 
bulge/BH growth. 


plenishment means widespread SF activity. Fig. 6 displays 
SF fraction versus redshift among bulge-dominated galax- 
ies in all five CANDELS fields. The sample's properties 
(redshift, Mstar, and SFR) are compiled/estimated by Yang 
et al. (2019). Here, we classify SF versus quiescent using 
the method in $2.5. The SF fraction values are calculated 
for galaxies more massive than 10! Mo, above which the 
bulge-dominated sample is complete up to z ~ 3 (Yang et al. 
2019). 

From Fig. 6, the SF fraction for bulge-dominated galax- 
les is generally low (e.g., ~ 0.2-0.3 at z ~ 1-2). The low 


SF fraction suggests that bulge-dominated galaxies do not 
have prevalent strong gas replenishment, which is required 
to maintain their SF activity. In comparison, we also plot the 
SF fraction for galaxies that are not bulge-dominated (e.g., 
disky or irregular) in Fig. 6. These galaxies tend to have high 
SF fractions (e.g., ~ 0.8-0.9 at z ~ 1-2), indicating preva- 
lent gas replenishment among them. 

We caution that the argument above only qualitatively sug- 
gests that bulge-dominated galaxies have weaker gas replen- 
ishment than non-bulge-dominated galaxies. It does not rule 
out intermittent gas accretion among bulge-dominated galax- 
ies. Intermittent gas accretion might lead to sporadic star for- 
mation, also consistent with the low but non-zero SF fraction 
among bulge-dominated galaxies (see Fig. 6). A quantitative 
assessment of intermittent gas accretion is beyond the scope 
of this work. 
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Figure 6. Fraction of star-forming galaxies (Mstar > 1019-2 Mo) 
as a function of redshift. The red data points represent bulge- 
dominated galaxies, and the blue ones are for the other galaxies 
that are not bulge-dominated. The error bars represent lo bino- 
mial uncertainties calculated with the "binom conf interval" func- 
tion of ASTROPY. The star-forming fraction is low (~ 0.2-0.3 at 
z œ 0.5-2) for bulge-dominated galaxies, indicating that gas re- 
plenishment is not common for them. 


4. SUMMARY AND FUTURE PROSPECTS 


In this work, we have presented the ALMA observations 
of 7 bulge-dominated star-forming galaxies at z = 0.5-2.5. 
Our main results are summarized below. 


* We have reduced the ALMA data and measured the 
CO(2-1) or CO(3-2) fluxes (see $2). We have detected 
the CO lines from 4 sources at > 3o levels, and we 
have estimated 3c upper limits for the other sources 
(82.3). From these results, we have inferred molecular 
gas masses (or upper limits) assuming aco = 4.3. By 
fitting the existing multiwavelength data with CIGALE, 
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we have estimated the stellar masses and star formation 
rates for our ALMA targets ($2.5). 


The gas-depletion timescales are at least 2 times 
shorter than the corresponding bulge/BH growth 
timescales for most sources, except for CAN- 
DELS 25998 (see $3.1). The median value of u is only 
< 0.26. If we assume aco = 0.8 (a typical value for 
compact galaxies) instead of aco = 4.3, all sources 
have u < 0.25 (i.e, Taep is > 4 times shorter than 
Tgrow) With a median « 0.048. The ALMA contin- 
uum measurement from Barro et al. (2016) also sug- 
gests u < 1 for two bulge-dominated sources (includ- 
ing CANDELS 25998). Therefore, we conclude that, 
without strong gas replenishment (supplying Mz; by 
a factor of = 4), the observed gas content of the SF 
bulges is generally insufficient to support significant 
bulge/BH growth. 


To assess gas replenishment, we have estimated the 
SF fraction for a mass-complete sample of CAN- 
DELS galaxies (see $3.2). The SF fraction for 
bulge-dominated galaxies is much lower than that for 
non-bulge-dominated galaxies (e.g., ~ 0.2-0.3 versus 
= 0.8-0.9 at z = 1-2). The low SF fraction of bulge- 
dominated galaxies indicates that gas replenishment is 
not a common process among them. We caution that 
our qualitative argument cannot rule out weak intermit- 
tent gas accretion among bulge-dominated galaxies. 


Our overall results indicate that the Mpy-Mstar rela- 
tion has already formed at the beginning of the bulge 
evolution phase (Scenario 2 in Fig. 1). The systems 
then maintain this relation until z = 0. In other words, 
the BHAR-bulge SFR relation has the role of main- 
taining the BH-bulge mass correlation, but not creating 
it. Therefore, it will be useful to study the BH-galaxy 
coevolution in the pre-bulge phase, which might re- 
veal the mysterious origin of the Mpy-Mstar relation. 
Such a study requires reliable techniques, probably 
machine learning trained by hydrodynamical simula- 
tions (e.g., Huertas-Company et al. 2018), to select 
pre-bulge samples from high-resolution images. 


Finally, we note that our sample size is limited (7 sources), 
and they are all below z — 2.5. If this small sample is some- 


how biased toward the late stage of star formation (82.2), 
then our main conclusion could be altered. Future (sub)mm 
observations of a much larger (Z; 100) bulge-dominated sam- 
ple over a wider parameter space (especially at z > 2.5) 
will naturally address this potential issue and further test 
our conclusion. High-redshift morphological classifications, 
which are necessary to select bulge-dominated sources, will 
be available in the near future with the advance of JWST ex- 
tragalactic surveys. ALMA or other (sub)mm facilities can 
perform follow-up observations of the JWST-selected tar- 
gets. 
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APPENDIX 


A. SPATIAL ANALYSIS 


For the 4 CO-detected sources ($2.3), we further perform 
spatial analysis for their lines. For each of the 4 sources, we 
divide the line into blue and red halves and make a line im- 


age for each half. The resulting line-image contours are dis- 
played in Fig. 7. From this figure, only CANDELS 24210 
appears to have (slightly) separated blue versus red con- 
tours. The angular distance between the blue and red peaks is 
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& 0.5". On the other hand, the positional uncertainty for the 
blue/red line emission is ~ 0.60(S/N) ! ~ 0.26", where 
0 ~ 1.5" is the synthesized beam FWHM (e.g., Ivison et al. 
2007). Therefore, the ~ 0.5" separation is marginally signif- 
icant at a ~ 2c level. 

The results above indicate that our current ALMA data are 
not able to spatially resolve the line-emitting regions in gen- 


eral. This is understandable considering the relatively large 
beam sizes compared to the H-band profiles (see Fig. 7). 
High-resolution ALMA runs are necessary to probe the spa- 
tial distribution of the CO emission. 
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